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Abstract 

Flexible boundary condition methods couple an isolated defect to bulk through the bulk lattice 
Green's function. The inversion of the force-constant matrix for the lattice Green's function requires 
Fourier techniques to project out the singular subspace, corresponding to uniform displacements 
and forces for the infinite lattice. Three different techniques — relative displacement, elastic Green's 
function, and discontinuity correction — have different computational complexity for a specified 
numerical error. We calculate the convergence rates for elastically isotropic and anisotropic cases 
and compare them to analytic results. Our results confirm that the discontinuity correction is the 
most computationally efficient method to compute the lattice Green's function. 

PACS numbers: 
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I. INTRODUCTION 



Atomic-scale simulation of isolated defects with a computationally tractable number of 
atoms requires careful choice of boundary conditions. Periodic or fixed boundary conditions 
introduce fictitious forces when relaxing the geometry of defects; reducing the error requires 
increasing the number of atoms. Flexible boundary condition methods avoid these errors 
by instead using harmonic lattice response for atoms away from the defect. In particu- 
lar, the bulk lattice Green's function (LGF) gives the short- and long-range displacements 
in response to a point or line force. Sinclair et al. [l| introduced flexible boundary condi- 
tions for studying defects such as cracks^, Q], dislocations , vacancies and free 
surfaces js] with classical potentials and isolated screw or edge dislocations with density- 
functional theory [9], Q, 11, 12]. Evaluation of the LGF in real space involves the inverse 



Fourier transform of a function with a singularity at the F-point {k = 0), which requires 
algorithmic approaches to evaluate numerically. Relative displacement method}^, [isl. elastic 
Green's function (EGF) correction [14] and discontinuity correctionjl^ are three techniques 
to numerically evaluating the bulk LGF. We compare all three methods to determine the 
most computationally efficient approach. Section [ITl defines the harmonic response func- 
tions, the relative displacement method, elastic Green's function correction and discontinu- 
ity correction, and the convergence evaluation methodology. Section IIIII follows with the 
convergence results and discussion. We find that the discontinuity correction has the fastest 
convergence rate over the relative displacement and elastic Green's function correction and 
verify our predicted convergence rates with a simple model and density functional theory 
results for Al. 

II. BACKGROUND 

The lattice Green's function G^{R — R!) relates the displacements u{R) of atom R to the 
internal forces f{R') on another atom R' of the lattice through 

u{R) = Y,G\R-I^)f{R')- (1) 

R' 
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Conversely, the forces on an atom can be expressed in terms of displacements through the 
force constant matrix Dj^R — R') by 

f{R) = -Y,DiR~li)u{R'). (2) 

R' 

Translational invariance of an infinite lattice makes a function of the relative positions 
of two atoms. Substituting eqn. ^ into eqn. ^ gives T^j^, G^iR - R')D{R') = -16{R), 
where 6 (jij is the Kronecker delta function. A constant shift in the atoms positions does 
not produce internal forces, giving the sum rule 'ZrIHR) = and making G^{R) the 
pseudoinverse of D_{R) in the subspace without uniform displacements or forces. Fourier 
transform of the lattice functions are defined as 

^ Jbz V'T^r 

R. 

for k in the Brillouin zone (BZ). The integral can be approximated by a discrete sum of Nk 
points as G^{R) = Y2k ^~^^^G^ {^) ■ In reciprocal space, the matrix inverse relation and 
the sum rule are G^{k)D{k) = 1 and D{0) = respectively. For a single atom crystal basis, 
D{k) expands as for small k as 

D{k) = - + ■ ■ ■] =^ E ■ ^)'^(^)' (3) 

R ' R 

due to the inversion symmetry of D.iR)- At the F-point, D(k) is of the order A;^, so G^(k) 
has a second order pole. The discrete inverse Fourier transform of G^{k) does not converge 
due to this singularity. 

Fig. [T] shows the relative displacement method, elastic Green's function correction 
and discontinuity correction which are used to avoid the singularity in LGF. Relative 
displacement. 



131 ] Rigid body translations leave the potential energy of the lattice un- 
changed, so it is enough to calculate only the relative displacements of atoms. Choos- 
ing an arbitrary atom as an undisplaced origin requires calculation of G^{R) - G^(0) = 
/ jSj^Q^i^y'"'^ - I j0fG^iky^-^ which reduces to 

G^{R) - G^(0) = j G^{k){cos{k-R) - l)dk, (4) 

due to the sum rule and inversion symmetry. For small fc, cos{k-R) — l is of the order k"^ which 
cancels out the second order pole in G^{k) leaving a /c-direction dependent discontinuity at 
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FIG. 1: Integrand of the inverse Fourier transform for LGF, relative displacement method, EGF 
correction and discontinuity correction at i? = (1, 1). G^(k) has a second order pole at the F- 
point. The relative displacement method avoids the pole by considering only the displacements 
relative to a fixed point. The EGF correction removes the second order pole by subtracting a 



at 



cutoff elastic Green's function. Removal of the pole creates a discontinuity independent of 
the F-point. The discontinuity correction removes the discontinuity created by EGF correction. 
The remaining part of the integrand is smooth in the entire Brillouin Zone. The bottom row shows 



the variation of the integrand as a function of 



when the origin is approached from different 



angles 9 = tan ^{ky/kx). A discontinuity at the F-point is created by the relative displacement 



but depends on the direction 



method and EGF correction where the integrand is independent of 
of approaching the origin. Notice the difference in vertical scale for the LGF — which has a pole — 
from relative displacement and EGF correction — which have a discontinuity — and the discontinuity 
correction — which approaches zero quadratically. 

the F-point. The discretized version of eqn. (jl]) is -^"^kicos^k ■ R) — l)G^{k). Elastic 
Green's function correction. Following the procedure and notations of [12] and using 



4 



— * 

eqn. ([3]), G (k) for small expands as 

= [eK^^\k) - k^A^^\k) + 0{k^)]-^ 

= k-^[A^^\k)]-^ + [A(2)(^)]-iaW(A;)[a(2)(^)]-i + Oik-'') 

= G''{k) + G^\k) + 0{k'), 

where k'^K'^'^^k) and k'^K^'^\k) are the second and fourth order terms in a small k expansion 
of D(k). The Fourier transform of the elastic Green's function G^(k) is the second order pole 
and G'^'^{k) is a /c-direction dependent discontinuity [15|]. The lattice Green's function can be 
separated into an elastic part G^{k) which should be inverse Fourier transformed analytically 
and the rest which no longer has a pole (it still has a discontinuity). The remaining part of 
the LGF can be inverse transformed numerically by J2k cos(A; ■ R){G^{k) — G^{k)fcut{k)) 
where fcut is a cutoff function that smoothly vanishes on the Brillouin zone edges. Removal of 
the second order pole by subtraction of a cutoff version of elastic Green's function is used in 



the semicontinuum method of Tewary IJ]. Discontinuity correction. To further improve 
convergence, the discontinuity correction treats the G'^'^{k) part analytically 15(]. In this case, 



the remaining portion of G^{k) given by Y.k cos(^ " R)iQ^{k) - {G^{k) + G'^''{k))f^ut{J^)) 
is smooth and can be integrated numerically more efficiently. 

We expect the convergence rate of the discontinuity correction method to be consistent 
with the results for integration of smooth periodic functions, while the convergence of relative 
displacement and elastic Green function correction methods should be dominated by the 
discontinuity. With N^i^ denoting the number of partitions in each direction, mid-point 



rule gives a A^^^^^ scale for convergence rate of such integrals in all dimensions 16|, [l7| . The 
number of k points Nk is Nf:^^ for dimensionality c? = 1, 2, 3; therefore, the convergence rate 
of the mid point rule scales as Nj^^^'^. However the EGF correction and relative displacement 
method should have a poorer convergence compared to discontinuity correction due to the 
discontinuity that they create at the F-point. Since the integrand is smooth elsewhere, we 
expect the error to be dominated by the area/volume around F-point and therefore be of 
the order of N^^ or N~^. 

We check the predictions of convergence for the three methods using (1) a simple cubic 
nearest neighbor model and (2) fee Al. First, as a simplified case we consider a square (cubic 
in 3D) elastically isotropic lattice with nearest neighbor interactions and lattice constant 



Oo = vr. We consider only one component of the full matrix: G^{kx,ky) = [sm'^{7ikx/2) + 
sin^(7rA;j^/2)]^^. The second order pole is given by the elastic Green's function G^{kx, ky) = 
— which is multiplied by a cutoff function to vanish smoothly at the BZ edges. The 



discontinuity correction is given by G {k^, ky) 
function. In three dimensions we have 



3 A: 



TTpr which is also multiplied by the cutoff 



G ikx-i ky, kz) 



G (kx, ky, kz) 
G {kx, ky, kz) 



[si\i^{Tikx/2) + sm^{'Kky/2) + sm^{'Kkj2)Y^ 
4 



TT 



3 k 



R) from 



121. The 



For our elastically anisotropic Al lattice, we obtain the force constant matrix D 
DFT using ultrasoft pseudopotentials with the generalized gradient approximation 
numerical integration over the BZ is done with a uniform mesh evaluating the integrand 
at the mid points. Even and odd values of A^div give meshes that include or avoid the F- 
point — what we call F and non-F centered meshes respectively. When applying the relative 
displacement method and EGF correction, the value of the integrand — which is discon- 
tinuous at = — is assigned zero at the F-point. We calculate the numerical error as a 
function of A'^^ and N^^^ to compare the efficiency of the three methods. 



III. RESULTS AND DISCUSSION 

Fig. [2] shows the convergence rates of relative displacement method, EGF correction and 
discontinuity correction in the square lattice case. The discontinuity correction and EGF 
correction scale as A"^^ and A^^^ respectively as expected. The value of R does not affect 
the power law scalings of the convergence. The prefactors on the other hand, are changed 
in the relative displacement method and are of the same order in EGF and discontinuity 
corrections. While the A^^^ convergence for relative displacement method obtained by a 
F centered mesh is in accordance with the analytical predictions, use of a non-F centered 
mesh produces a convergence faster than expected for this method. This is an artifact of 
the isotropy of the EGF. 

The integrand in the relative displacement method is I{k) = (cos(/c ■ R) — l)G^{k). Near 
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FIG. 2: Convergence rate with number of k-points of the relative displacement method, EGF 
correction and discontinuity correction in a 2D square lattice. We expect N^^"^ convergence for 
discontinuity correction and poorer N^"^ convergence for EGF correction and relative displacement 
method. Using a non-F centered mesh (left) causes an unusually fast convergence for the relative 
displacement method in elastically isotropic materials. The exponents in the power law scalings 
are not affected by the value of i?, while prefactors are changed in relative displacement method 
and are of the same order in EGF and discontinuity corrections. 

the F-point, G^{k) matches G^{k) and the leading term in the integrand is 
For an isotropic EGF, G^{k) is constant, so 



r , . rko/'2 rko/2 qE ,2 

li2ik)d^k= / {-—R''cos\9))dKdky = -^G''R\ (5) 

Jk^ J-ko/2 J-ko/2 2 4 



where ^ is the angle between vectors k and R. The value of the integral over a square 
ko X ko region around /c = 0, for small k is 

r.fco/2 /•ko/2 nE 72 
-fco/2 J-ko/2 

The midpoint rule integration of the same region with a non-F centered mesh uses the k 
points h = iko/2,ko/2), h = {-ko/2, ko/2), h = {-ko/2, -ko/2) and h = {ko/2, -ko/2) 
each contributing area k^/A. The angle between each ki and R, are 61, 62 = 61 + 11/2, 
6^ = 61 + n and 64 = 61 -\- 37r/2. Therefore, the numerical approximation for the integral 
around the F-point is 



/r2 

f _ 

liso - ^ 



€\h) + lZ\k2) + €\h) + €\h) 

k^ ' k^ 

^G^R' [cos^ e, + sin^ 9, + cos^ 9^ + sin^ 9^] = -'^G^R\ 
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which is equal to the exact value of the integral around F-point given by eqn. (I5j). To avoid 
the effect of the discontinuity at the origin using a F centered mesh, the F-point contribution 
to the integral is considered zero while its actual value is given by eqn. ([5]). This is the source 
of the dominant error in relative displacement method on a F centered mesh which according 
to eqn. ([5]) accounts for the R dependence of the error. The R dependence of the error is 
verified by comparing the ratio of prefactors of the relative displacement convergence laws 
for different R values and the corresponding i?^ in Fig. [2] which are both approximately 
27. On the other hand, the non-F centered mesh automatically gives the exact value of 
the integral around the origin based on eqn. ([6]) and thus produces a faster convergence 
limited only by the convergence of smooth periodic functions. Note that if G^{k) depends 
on k — which is the case for anisotropic elastic response — the numerical approximation of 
the integral around F will not be equal to its exact value. Therefore, in general anisotropic 
problems the non-F mesh is not special. 

Fig. [3] shows that the 3D results follow the same trend as the 2D ones in accordance with 
the expected values. Both F centered and non-F centered meshes give 

AT-^/-^ (rf = 3) and A^-i 

scale for the convergence rate of discontinuity correction and EGF correction respectively. 
Similar to the trend observed in 2D case the F centered mesh produces the expected N^^ 
scale for the convergence of relative displacement method and the non-F centered mesh 
produces faster convergence due to the isotropy of the elastic Green's function. 
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FIG. 3: Convergence rate with number of k-points of the relative displacement method, elastic GF 
correction and discontinuity correction in a 3D cubic lattice. The error for discontinuity correction 
method scales as N^^ '^^'^ where the dimension d is equal to three here. Note that using a non-F 
centered mesh creates a faster convergence for the relative displacement method as observed in the 
2D case. 




Nj^(# of k points) N^(# of k points) 



FIG. 4: Convergence rate with number of k-points of the relative displacement method, elastic 
Green's function correction and discontinuity correction in computation of the Gn component of a 
2D LGF in Al. The Al lattice constant ag is 4.04A. The convergence trend of LGF calculations in a 
FCC lattice is the same as the one observed in the simplified problem which is also consistent with 
the expected values. Note that use of the non-F centered mesh does not cause a fast convergence 
for relative displacement method due to the anisotropy of the elastic Green's function. 
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TABLE I: Effect of dimension on the convergence rate with number of k-points and number of 
divisions for the relative displacement method, EGF correction and discontinuity correction. A'div 
is proportional to , the inverse grid spacing and = -^div- "^^^ discontinuity correction scales 
as -/V^^ (or A'^^ ^^^) while the EGF correction and relative displacement method scale as (or 



2D 3D 



Power law scaling of error with 




A^div 


Afc 


Adiv 


Disc correction 


-2 


-4 


-4/3 


-4 


EGF correction 


-1 


-2 


-1 


-3 


Rel. displacement 


-1 


-2 


-1 


-3 



Fig. m shows that the convergence trends are not changed for an anisotropic long range 
interaction — fee Al — except for relative displacement method. The convergence rates of the 
three methods for the two dimensional LGF are shown. This is a relevant case that occurs 
in modeling dislocations. The lattice is periodic in the threading direction [110] which is 
appropriate for studying screw dislocations. With a non-F centered mesh, the anisotropy 
of the elastic Green's function eliminates the fast convergence of the relative displacement 
method. The convergence trends of the three methods show that these trends are not specific 
to the simplifying assumptions of isotropy or short range interactions and therefore can be 
trusted in realistic calculations. 

Table [T] summarizes the convergence results for the three methods. The expected con- 
vergence rate for a numerical integral of a smooth periodic function evaluated by mid-point 
rule is Ajr^. When expressed in terms of the number of A;-points used in evaluating the 
integral Nk, the convergence rate would be proportional to A^^ Since the discontinuity 
correction leaves a smooth periodic part of the integrand, it follows the above convergence 
rate. The EGF correction and relative displacement method also converge with the scale 
of A^^^ or A"^^. Therefore the discontinuity correction method has the fastest convergence 
rate. The convergence rates imply that a certain amount of error is achieved with less A^^ 
by discontinuity correction method compared to EGF correction or relative displacement 
method which means that the discontinuity correction requires the least computational ef- 
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fort. Although the EGF correction and relative displacement method require comparable 
computational effort, the R dependence of the pref actors suggests that the relative displace- 
ment method takes even more fc-points than the EGF correction. Also note that there 
is a trade-off between less computational effort and more complex algorithms. EGF and 
discontinuity corrections calculate the elastic Green's function and discontinuity correction 
parts of the LGF analytically while relative displacement method does not require additional 
analytic evaluations. 

IV. CONCLUSION 

We find the most efficient method to compute the lattice Green's function to be the dis- 
continuity correction. The relative displacement method, elastic Green's function correction 
and discontinuity correction have all been used in different calculations; we applied and 
compared the three methods to calculate LGF for a square nearest neighbor lattice and long 
range fee Al. The convergence trends match the analytical values with an unusual exception 
for lattices with isotropic elastic Green's function. It is shown that the discontinuity correc- 
tion improves the convergence rate to quadratic convergence for 2D calculations compared 
to linear convergence for the relative displacement and elastic Green's function correction. 
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